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Instantaneous Frequency and 
T ime- Frequency Distributions 

Boualem Boashash and Graeme Jones 



Keywords: instantaneous frequency, group delay, bandwidth- 

duration product, analytic signal, Hilbert transform, quadrature 
signal, digital instantaneous frequency 



1 Introduction 

This chapter presents the fundamental principles on which the notions 
of instantaneous frequency (IF), group delay (GD) and analytic signal 
are based, in the light of their relationship with time-frequency 
distributions (TFDs). An understanding of these notions is necessary 
in order to efficiently use them to derive methods for time-frequency 
analysis of non-stationary signals. 

As indicated in chapter 1, many natural and man-made signals 
have spectral characteristics which vary in time. The chirp (linear 
FM) signal is such an example: it consists of a sine wave whose 
frequency increases (or decreases) linearly with time. Such a signal is 
used in practical situations such as seismic surveying, communications, 
radar, sonar. It also occurs naturally in such places as the echolocation 
systems of bats. 

For the chirp and other time-varying signals the instantaneous 
frequency (IF) is an important descriptor. Qualitatively, this is a 
parameter which corresponds to the frequency of a sine wave which 
locally (in time) fits the signal under analysis. It has meaning only for 
monocomponent signals, that is, where there is only one frequency or 
a narrow range of frequencies varying as a function of time. In other 
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cases, the notion of a single valued IF becomes meaningless, and some 
form of appropriate break-down into sub-components is necessary. 



2 Basic concepts 

2.1 Evolution of the term 'instantaneous 
frequency’ 

This section presents a review of some early work concerned with 
defining and interpreting the IF. Carson and Fry in 1937 [1] generalised 
the notion of frequency by allowing it to vary as a function of 
time. They argued that it should be termed the IF, since it was a 
generalisation of the definition for constant frequency i.e. the rate of 

change of the phase angle at time t. 

In 1946 Van der Pol [2], in the context of communication theory, 
defined the IF by considering a signal of the following form: 

s(t) = a 0 cos(27 rft + </>o) (1) 

where ao is a constant, / is the frequency of the oscillation, and <j>$ is 
a phase constant. 

In this expression, Van der Pol makes ao and </> o vary and defines 
the amplitude and phase modulations as: 

a(t) = ao 1 T kh(t) (2) 

a Y~i ri 

0(t) = </>o[l J r m 9(t)} , ( 3 ) 

where k and m are constants, and h(t) and m(t) are the time- varying 
modulating waveforms. One could then express the frequency as: 



/ = /o[l + rng(t)\ ( 4 ) 

Equation (4), however, leads to physical inconsistences as there is 
time- varying information in the modulated phase. This may be seen 
by the substitution of / in (4) back into equation (1). 

Van der Pol argued that to properly obtain an expression for the 
IF of a phase modulated signal, equation (1) should be re-written as: 



s(t) = a 0 cos [ 27 vfi(t)dt + fa 

Jo 



To further explain the concept, Van der Pol used an analogy with 
harmonic motion. For a vector rotating at constant angular velocity 
the frequency is: 



/ 



1 d(j){t) 
27r dt 



(6) 
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where <f>(t) = 2 tt ft is the phase, and may be considered to be the angle 
the rotating vector forms with the horizontal reference line. He argued 
that for a varying angular frequency (represented by a vector rotating 
with a variable angular velocity), one could still use the definition 
given in (6) and obtain an (instantaneous) frequency of rotation: 



MO 



1 dcj)(t) 
2tt dt 



This vector may be visualized as the complex signal (see figure 2.1): 

= ( 8 ) 

which is the sum of the real signal in (5) and the pure imaginary 
quadrature signal: 



ja(t) = ja 0 sin / 2tt fi(t)dt + <j> Q 

Jo 






a r 



0(t) 



Figure 2.1. A vector (complex signal) rotating with variable angular velocity (instantaneous 
frequency). 

Thus the phase of the resultant complex signal is used to define 
the IF of (7) which is consistent with the generally accepted definition. 

Gabor [3] proposed the Hilbert transform as a means of obtaining 
a unique complex signal from a real one — with the imaginary part 
in quadrature to the real part (equation 9). The Hilbert transform is 
defined by: 



7i[s(t)) = p.v. 



+°° s(t ~ t) 



( 10 ) 
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where p.v. denotes the Cauchy principal value of the integral. Ville 
4] used this transform to generate his complex analytic signal. With 
the IF as defined in (7), he showed that the average frequency (using 
Gabor’s average measures [3]) in a signal’s spectrum was equal to the 
time average of the IF [4]: 



where 



< / >=< fi > 




( 11 ) 

(12) 



and 




The analytic signal is z(t) = s(t) -\-H [$(£)], and Z(f) is the Fourier 
transform of z(t). 

The interpretation of the IF has been a controversial issue. For 
example, Shekel [6] argued that the IF defined by: 



Mt) = 



1 dcf>(t) 
27 r dt 



(14) 



where z(t) = a(l)e-^W is not a unique function of time, since,, for 
example, any AM wave (written in complex form) may be expressed 
as either m(t)e j2nft or m 0 e^W, where m 0 is a constant. A unique 
complex representation is obtained by using the HT, as Gabor and 
Ville have stated [3] [4], but whether or not it corresponds to any 
physical reality is another question, which is discussed in a later 
section. 

Mandel [7] also challenged any physical interpretation of the IF by 
providing examples where the IF of a signal does not correspond to 
any of its Fourier spectral components. In addition, he showed that 
for moments higher than the first, the Fourier frequency and the IF 
average measures do not coincide, as Ville had noted for the second 
moment [4]. Several other authors have also contributed to the study 
of IF in more recent work [8] [9] [10] [11 . 



2.2 The Hilbert transform and the analytic 
signal 

Gabor first used the Hilbert transform (HT) to generate his complex 
signal according to: 
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z(t) = s(t) +jy(t) (15) 

where y(t) represents the HT of s(t), as defined in (10). This signal 
is referred to as the analytic signal. It does not always correspond, 
however, to a signal plus its imaginary quadrature component. As an 
illustration, suppose we have a signal of the form: 

s(t) = Hr(t) c ° s ( 27 t fot + </> q ) ( 16 ) 

where Ur(t) denotes a box function of unit amplitude, duration T and 
symmetric about t — 0. 

By inspection, it is obvious that the signal in phase quadrature is: 

q(t) = n r (t) sin(27 rf 0 t + (j>o) (17) 

Applying the HT to the signal of (16) does not exactly yield the 
signal in phase quadrature (equation 17). This is because the HT 
operation preserves the positive frequency domain of the spectrum 
and inverts the sign of the spectrum in the negative frequency domain. 
Thus in this case it does not simply transform the cosine term into a 
sine term. 

It is therefore essential to know when the HT generates the signal’s 
quadrature component exactly. This can be determined by considering 
a signal of the form a(t) cos <£(£), and defining conditions for the 
following equation to be verified: 

a(t) cos (j>{t) + jli[a(t) cos 4>{t)] = a(t)e 

— a(t) (cos <j)(t) + j sin </>(£)) 

Equation (18) is satisfied if one of the following two conditions 
holds: 

1 The spectrum of a(£), A(f ), lies entirely in the region I/I < fo 
and jF{cos only exists out of this region [12] [13]; and 

2 A(f) is non-zero only for / > — /i , and jF{cos />(£)} is non-zero 
only for / > /2 ; for > /i > 0 [4] [14] [35] where T denotes 
the Fourier transform operation. 

These two conditions are obtained by examining the result of the 
convolution of A(f) and .F{cos </>(£)}, as T{a(t) cos </>(£)} = A(f) * 
J 7 {cos <p(t)} . It should also be stated that condition 2 above is not a 
practical constraint as a(t) is defined to be an amplitude envelope and 
thus it must be low frequency, centred at the DC point. 

The HT-generated analytic signal is meaningful and representative 
of the physical signal only under some conditions. As explained above, 
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if we have a signal with an imposed modulation of the form a(t) cos </>(£) 
where physical meaning is attached to a(t) and. (f>ft), and if the spectra 
of aft) and cos <f>(t ) are not separated in frequency, then the HT will not 
be able to distinguish the amplitude and phase functions. Although 
the analytic signal will be of the form: 

(19) 

and will be unique, a z (t) and </> z (t) may have no practical meaning. 
The more closely a signal approaches a narrowband condition, the 
better is the analytic signal approximation to the quadrature signal. 



2.3 The bandwidth-duration(BT) product 

The previous work examining the HT, analytic signal and the IF has 
not considered practical aspects of signals. In practice, signals have 
imposed constraints such as causality, finite energy, finite duration and 
stability. Often, it is also assumed that they possess a finite spectral 
bandwidth (within which most of the signal energy is contained). 
Although a signal with both a finite duration T and a finite bandwidth 
B is prohibited by the Fourier transform, signals which approach this 
double limitation form an important class. A review of such concepts 
and measures for time-frequency bounds are given in 15]. 

When signals have a BT product sufficiently high, the approxi- 
mation error involved in assuming band and time limited functions is 
very small. These signals are often referred to as ‘asymptotic’ [16], 
the ‘asymptotic’ behaviour being measured by the parameter BT. 

With these notions in mind, the results involving the HT and 
the analytic signal from the previous section may now be put in a 
more general and quantitative framework, which is summarized in the 
following theorems: 



Theorem 1 

If the signal, s(t) = aft) cos (f(t), has a monotonic and positive 
frequency law, then the quantity sft ) + jhi[sft)\ approaches z(t) — 

aft)e MW asymptotically as BT 00 . 

The proof follows work by Bedrosian [12] and Vakman [17] and is 
based on the following relation holding when BT 00 . 




H[s(t)} = H[a(t) cos <t>(t) 



= a(t)H [ cos <j>{t) 



= aft) sin <fr(t) 




It should be noted that the theorem states that the IF must be 
positive for the duration of the signal, and it must be monotonic in 
time [18]. 

Theorem 2 

For a signal of the form, sft) = aft) cos(27r/ot + </>(t)) (fo i s the 
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constant or tonal frequency, not present in <f>(t) ), with an arbitrary 
BT value, the analytic signal derived using a Hilbert transform 
asymptotically approaches the signal plus its imaginary quadrature 
part, i.e., a(£)e J *( 2 ^° t+< ^), as fo — * 00 . 

The proof follows work by Nuttall [14], where he shows that the 
difference in energy between the HT generated analytic signal and the 
signal plus its imaginary quadrature part is the energy in the spectrum 
of: 



So(f) = r{a(t)e m) } (21) 

for / < -fo. Thus the representation becomes exact as /o — > 00 . 

Theorem 2 indicates that the approximation error in forming the 
complex signal with the HT (compared to the signal plus its imaginary 
quadrature component) and the BT product are not directly related. 
In effect, theorem 1 is a special case of the more general theorem 2. 
A good approximation to the a(t)e^ 27r ^° t4 '^ ) ^^ term (and thus the 
IF) can be achieved even for a low BT product signal, provided its 
centre frequency is large. Two examples are used here to illustrate 
the theorems. 



Example 1: Consider two signals of the form, 



sft) = n T (t)cos 27 r 





Signal 1: / 0 = 50Hz, a = 1714.3, BT = 2.1 
Signal 2: f 0 = 5QHz, a = 93.7, BT = 38.1 



( 22 ) 



Figure 2.2 shows the (analytic) spectrum of signal 1 (created using 
the HT). This may be compared to figure 2.3 which is the spectrum 
of signal 1 plus its phase quadrature component. Figures 2.4 and 2.5 
show the same results for signal 2. Examination of the figures provides 
a simple illustration of theorem 1 — the analytic signal approaches 
the real plus quadrature signal as the BT becomes large (this may be 
seen for the spectra of signal 2). 



Example 2: Consider two signals, once again, of the general form of 

(22): 



Signal 1: fo — 11Hz, a — 333.3, BT = 1.2 
Signal 2: fo = 51Hz, a = 333.3, BT = 1.2 



Amplitude 
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Analytic signal spectrum (BT = 2.1) 




Figure 2.2. The analytic signal spectrum for signal 1 of example 1. 
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In a similar manner to the first example, figures 2. 6 (2. 8) and 
2. 7(2. 9) give the analytic spectrum and the real plus quadrature 
spectrum respectively for signal 1 (2). From the figures it is obvious 
that the analytic signal approaches the real plus quadrature complex 
signal as /o becomes larger, which illustrates theorem 2. 



2.4 Instantaneous frequency and group-delay 

Let s(t) be a continuous time real signal with z(t) its corresponding 
analytic signal: 



TgU) = 



1 cMJ) 

27 r df 




■ l n many applications it is common to use the GD, r p (/), of the 
signal, to characterise the time-frequency law of the signal. It is 
therefore desirable to relate the two quantities of IF and GD. 

In order to achieve this, the Fourier transform of the signal is 
examined. For signals of the form, a{t)e^\ with a large BT product 
and a monotonic IF law, the Fourier transform may be approximated 
by application of the stationary phase principle as follows [17]: 



z(t) = a(t)e?+® (23) 

The IF is simply obtained by differentiation of the phase and 
scaling by z(t) has a complex spectrum which may be expressed 
as: 



Z(f) = A(f)e m (24) 

where a(t) and A(f) are positive functions. Another quantity of 
interest for the signal is the group delay (GD) defined by: 



J- {a{t)e^^} « 



2ix 



a(A)e^“ 2 ^ A+ ^ A ) ±7r / 4 l 



where A is the stationary phase point, that is, the solution of : 




d 

dt 



(-j2irft + <t>(t)) = 0 




It is simple to show, using equation (26), that the IF and GD are 
approximately inverses of each other (i.e., fi(t) = r~ 1 (/)) [18] for these 



Amplitude 
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Analytic signal spectrum (BT = 38.1) 
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Figure 2.4. The analytic signal spectrum for signal 2 of example 1. 

large BT signals. A further discussion of the relationship between IF 
and GD (including analytical examples) may be found in [15]. 

Signal concentration and elementary cells 

A related concept utilizing IF and GD laws is that of ‘ elementary cells 
which define the signal concentration about these laws. One way to 
define the dimensions of these cells is to introduce the notion of “re- 
laxation time”, T r , and “dynamic bandwidth, Bd , as [8]. 



T r (t) 



dm 1/2 

dt 



(28) 



(29) 



These measures represent boundaries in time and frequency within 
which most of the signal energy is contained. They provide quantita- 
tive support to the intuitive notion that the energy of an asymptotic 
signal is concentrated in a finite domain in the general time-frequency 

plane. 



Amplitude 
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Real + Quadrature spectrum (BT = 38.1) 




Frequency (Hz) 

Figure 2.5. The real signal plus quadrature signal spectrum for signal 2 of example 1. 



3 The IF and time- frequency distributions 



The following section examines the relationship between the IF of a 
monocomponent signal and its time-frequency representation. 

While working on the concepts of the IF and the analytic 
signal, and borrowing previous results from quantum mechanics, Ville 
formulated a distribution of a signal in time and frequency referred to 
as the Wigner-Ville Distribution (WVD) [4]: 



W(tJ) 



Too 



z(t + r/2)z*(t — r/2)e j27r ^ T dr 



(30) 



It was shown in [4] that the first moment of the WVD with respect 
to frequency yields the IF: 



m 



Too 



Too 



’ W(t, f)df 



(31) 



Ideally, one would expect from a time-frequency representation, of 




Amplitude 
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Frequency (Hz) 

Figure 2.6. The analytic signal spectrum for signal 1 of example 2. 



a signal with a large BT product, that the energy would concentrate 
about the IF, with a spread somehow related to the envelope of 
the signal. Accordingly, for monocomponent signals of the form of 
a(t)e it would be intuitively satisfying to generate a TFD that 

could be expressed as: 



p{tj) = A(tJ) & 6(f - m) (32) 

where A[t : j ) is the time-frequency representation of the amplitude 
function, a(£). Thus the amplitude and phase would be separable, 
providing an easily interpretable distribution: the distribution is 

centred around the varying IF about which the amplitude information 
is spread in time and frequency. 

Unfortunately, no TFE) exists that can be expressed in the form 
given in (32). For example the WVD (for a signal with a large BT 
product) yields a distribution of the general form [22]: 





W(t, f ) = kAi 



(33) 



Amplitude 
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Real + quadrature spectrum (BT = 1 .2, f 0 = 1 1 Hz) 




Frequency (Hz) 

Figure 2.7. The real signal plus quadrature signal spectrum for signal 1 of example 2. 

/ •Too 

W(t 1 f)df = 

-00 

a 2 (t). 

Only for the case where the third order derivative of the phase 
function is zero (linear FM signals), does the WVD produce a 
distribution of the form of (32). In other cases, where the phase is 
more complicated, the WVD is distorted by the Fourier transforms of 
higher order phase derivative terms. 

Many popular TFDs may be calculated by a single generalised 
formula, derived by Cohen [23]: 

/ Too c-foo c+00 

/ / ef^-Vg&T) 

-OO J—OO J —OO 

(34) 

z(u + t/2)z*(u — r / 2)e~~ J 2?r ^ T dududr 

where g(v,r) is the kernel function which defines the particular TFD 
chosen. The IF will be given by the first frequency moment of any 

TFD satisfying (34), provided the following conditions hold [15] [23]: 




Amplitude 
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Real + quadrature spectrum (Bt = 1 .2, f 0 = 51 Hz) 




Frequency (Hz) 

Figure 2.8. The analytic signal spectrum for signal 2 of example 2. 



/ x , dq(u.r) 

g(v, 0) = 1 and — = 0 

V ' dr r=0 

These conditions are in addition to those that must be satisfied by 
g(v,r) to be a member of Cohen’s class [24]. 

Many distributions yield the IF by correct first moment calculation 
as discussed above, but this is often computationally expensive and 
adversely affected by noise. It would be an advantageous feature of 
a distribution to allow estimation of the IF by peak detection. With 
signals whose phase functions are quadratic, the WVDs will be of 
the form of (32) and thus IF estimation can be achieved by peak 
detection. Such a method has been shown to be optimal for high to 
modest signal-to-noise ratios (SNRs) [24]. 

The magnitude squared short-time Fourier transform (STFT) has 
for a long time been the standard means of analysing non-stationary 
signals. It may be expressed as [15]: 

z(r)h(t 

where h(t) is the analysing window. 



Sz(tJ) I 2 



2 

r)e- j2llft dT (36) 




Amplitude 
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Analytic signal spectrum (BT = 1 .2, f 0 = 51 Hz) 




Frequency (Hz) 



Figure 2.9. The real signal plus quadrature signal spectrum for signal 2 of example 2. 



As the STFT is a straightforward extension of the stationary 
Fourier transform, the representation of a signal and its IF law is 
dependent on the length of the analysing window. It may be shown 
that for a signal whose IF law varies by an amount 6fi during the 
time interval 8t , the IF can be tracked (by peak detection) to within 
an accuracy A / provided that [18]: 

< (A/) 2 (37) 

where A / would be the spectral resolution width of the window, h(t). 
For a signal with a linear IF law, this reduces to a simple expression 
for the optimum window width [26] [27]: 



6fi 

8t 




This is the largest length window which will ensure that the IF 
is resolved as an unbiased peak in the STFT. The window must be 
small enough so that the IF does not vary too much (i.e., the IF is 
appioximately stationary during the analysis time) and large enough 

to provide the best possible spectral resolution. 
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3.1 The concept of a multicomponent signal 

The signals discussed so far are monocomponent. If the signal is 
composed of several components, the IF as previously defined does 
not always have immediate physical significance. Consider the case of 
a signal: 



s(t) = ai(t)cos0i(t) -F a 2 (t) cos fait) H (39) 

Assuming that each term of the signal s(t) has a large BT product, 
application of the HT to produce the analytic signal approximately 
gives: 



z m (t) — a z (t)e^ z ^ = ai(i)e^ 1 ^ F a 2 (t)e^ 2 ^ F • • 
The IF of this signal is: 



(40) 



n 



n 



l 



m) 



E^wFoo + E pij(fiid) + Fd)) 

i-l 



n 



n 



(41) 



E«*F + E «*(*) 

ij=l[i^j] 



i — 1 



where qij(t ) — ai(t)aj(t) cos(^(t) — <fj{t)) and the derivatives of the 
amplitude terms are assumed negligible. 

Little or no meaning can be attached to this IF value, which is 
oscillatory and highly non-linear, and it is in such situations that the 
concept of a multicomponent signal (and indeed a monocomponent 
signal) becomes attractive. Thus, this tentative intuitive definition of 
a multicomponent signal is proposed: 

Definition 

A finite energy asymptotic signal, s(t), is referred to as a multicom- 
ponent signal if there exists a finite number, N , of monocomponent 
signals, S{(t) ( i = 1 ,N), such that the relation s(t) = Y^i=i s i(t) 

holds for all values oft for which s(t) is defined, and this decomposi- 
tion is meaningful. 



Naturally, this definition must be applied with care as there can 
be no unique decomposition into individual monocomponent signals 
without a priori knowledge of the signal. 

A consideration of mono and multicomponent signals, in conjunc- 
tion with TFDs, leads to further interesting results and observations. 
The concept of a multicomponent signal is both a local and a global 
phenomenon. As figure 2.10 illustrates, the signal is multicomponent 
everywhere except where the two signals cross, and at that moment 
the signal is truly monocomponent. The concept of a multicomponent 
signal does have a global significance as the real signal must be known 



for all time to construct the analytic signal and the distribution. Ad- 
ditionally, the individual signals in a multicomponent signal do not 
have unique time-frequency laws unless the signal is always multicom- 
ponent (i.e., the signals’ laws do not cross). As an example consider 
the signals displayed in figure 2.10 — two possible decompositions into 
monocomponent signals are shown. 



t t 




Figure 2.10. Alternate ways of separating two frequency laws (decomposition of 
multicomponent signals is not unique unless the frequency laws of the components do not 
cross at any point). 



4 IF estimation via TFDs 

A number of approaches may be adopted for estimation of the IF of a 
non-stationary process. We present here the methods which use TFDs 
directly. 



4.1 IF estimation using the first moment of a 
TFD 



An approach to the estimation of IF, which uses moments of TFDs, 
and in particular, the WVD, was proposed in [30]. Taking the first 
moment of the general TFD expressed in equation (34) yields [15] : 
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m l p {t) 



where G(t, r) is the Fourier transform of the kernel function, g{v, r) 
of (34). 

This equation shows that the first moment of a TFD is, in general, 
a smeared version of the IF. If the signal magnitude is constant, then: 

m l p {t) = fi(t) (*) G(t, 0) (43) 

Thus in general the operation of estimating the IF by use of the 
first frequency moment of a TFD results in a smoothed version of the 
IF: 

fm = fid) * h d ) 

where h(t) — G(t, 0). 

For the case of the WVD, G(t, 0) — 6(t) and there is no smoothing, 
whereas for an STFT, h(t) becomes the original analysing window [31] 
and the result is a smoothed estimate of the IF, which when noise is 
present has a low variance but is characteristically biased [15] [39]. 
An examination of the STFT and its representation of IF, including 
results for choosing optimal windows, are given in chapter 8 of this 
book by Harris. 

Due to the high computational complexity, IF estimation utilizing 
the first moment of a TFD is not always a preferred method. 
The particular advantage of this technique, however, lies in the 
preprocessing that may be performed in the time-frequency plane (see 
section 5.1) — noise effects may be reduced as well as possibly allowing 
IF estimation of multicomponent signals. 



(44) 




4.2 Peak detection from TFDs for IF estimation 



Optimal estimation of the frequency of a sinusoidal process in 
Gaussian noise is achieved by peak detection from the magnitude 
squared STFT [33]. Where the signal is a non-stationary process, 



it would seem reasonable to employ a similar method to estimate 
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the instantaneous frequency (IF). The short-time Fourier transform 
(STFT), however, has poor energy concentration properties for rapidly 
time- varying signals [15]. IF estimation using the peak of the WVD 
may be used. The performance of this estimator, however, significantly 
degrades at low SNR [34]. The cross WVD (XWVD), which achieves 
high energy concentration as well as excellent noise preformance, was 
proposed for obtaining an IF estimate [18]. The XWVD between a 
reference signal z s (t) and an observed signal z r (t) is: 

/ °Q nr 

z s (t + T/2)z r (t-r/2)e-^f T dr (45) 

-OO 

The following recursive IF estimation procedure is then: 

1 Obtain an initial estimate of the IF [15] and form a unit 
amplitude signal that has this IF estimate as its frequency law; 

2 Using this as reference form a XWVD and extract the peak (of 
the squared XWVD) as the new IF estimate; and 

3 Repeat the procedure from step 1. 

Each time a new XWVD is estimated, the signal energy is more 
concentrated, so that there is a greater probability that its peak will 
be correctly estimated from a background of noise. For chirps at 
high SNR, the performance of the XWVD scheme will approach the 
performance of the STFT for estimating a corresponding sinusoidal 
frequency. As the SNR increases, the estimate will still exhibit very 
low variance compared with spectrogram estimates, provided the SNR 
and the spectral variation obey certain criteria. An error analysis and 
convergence performance are provided in [18]. As an example consider 
a linear FM signal in 0 dB white noise. The initial estimate of the 
IF (using peak detection from an STFT) is shown in figure 2.11. The 
estimate after one iteration of the above-described method is displayed 
in figure 2.12 (The true IF is shown in figure 2.13). It may be seen that 
the updated estimate is a significant improvement over the original. 
The algorithm typically converges in a few iterations. 



5 Instantaneous frequency estimation for 
discrete time signals 

This section provides a brief introduction to the problem of IF 
estimation for discrete time signals. A comprehensive study can be 
found in [39] and [42]. 
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Figure 2.12. The IF estimate of the reference signal obtained from the XWVD scheme after one iteration. 
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5.1 Discrete IF measurement s based on finite 
differencing 

The expression for the IF of a continuous time complex signal, z(t), 
may be re-written as: 

fid) = Kin ( 4>{t + St) - 4>(t - St)) (46) 

This may be translated directly to the discrete time case which 
provides an estimate of discrete instantaneous frequency (DIF) — 
where for example a central finite difference operation replaces the 
differentiation [19, part II]: 

fi(n) = ~(</>(n + 1) - 4>(n - 1)) (47) 

47T 

It is assumed here that the discrete phase, </>(n), has been 
unwrapped [39]. The discrete time analytic signal, z(n) J associated 
with the real discrete time signal, s(n), is given by: 

z(n) — s(n) + jH[s(n)} (48) 

where H[ ] is the discrete time Hilbert transform [20] which may be 
defined by: 

ff[ s („)] = V (modd) (49) 

Z— / m7r 

m =— oo 

Other estimates based on finite differencing have been proposed 
[21]. For a discrete time signal of the following form: 

x(n) = A cos <f){n ) (50) 

Griffith estimates the DIF as: 

fi{n) = ^(<A( n ) i)) (51) 

which is convenient for signals with modulation laws expressable as: 

M n ) - fo + rn(n) (52) 

Claasen and Mecklenbrauker indicate in [19, part II] that the 
central finite difference DIF of (47) is more appropriate than the 
backward difference estimator (equation 51) or the forward difference 
(4>(n 4- 1) - ^(n)) since the latter are asymmetrical. 

The central finite difference does not always produce unbiased 
estimates of the DIF. For IF laws with a non-linear variation, the 
estimate will be biased. An unbiased estimate for the DIF which may 
be expressed as a generalised finite difference estimator is discussed in 
[39] and [42]. 
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5.2 Zero-crossing measures 

A traditional and very simple method for estimating the DIF is 
the zero-crossing detector, which detects the discrete time distance 
between adjacent zero crossings, and then makes a frequency estima e 
assuming a stationary sinusoid: 



where T z is the zero crossing interval. .. 

This estimator is straightforward to implement, is naturally 
suboptimal and will produce biased estimates for non-stationary IF 
laws. Various forms of averaging may be used to reduce the estimator s 



variance 



[411 



5.3 Short-time estimates 

A number of estimation techniques have been proposed based on the 
assumption that the DIF is stationary over the local analysis time. 
As a straightforward extension of the periodogram, the magm u e 
sauared STFT provides an estimate of the DIF by peak detection. 
This estimate becomes maximum likelihood in high SNR for a linear 
DIF, providing the analysis window has been appropriately selected 

[391 . 

The resolution of such estimates may be improved by the use 
of parametric techniques. Griffiths [21] used a least mean square 
(LMS) linear prediction based spectral estimate. The prediction 
filter coefficients are recursively updated to reduce the computations 
required and preserve the local character of the estimates. The method 
is unable to resolve fast changing DIFs, but as a result it has good noise 
rejection properties for slowly varying laws. Other techniques, based 
on recursive least squares (RLS) approaches to the linear prediction 
problem, and time-varying autoregressive models, have been suggested 

[28]. 



5.4 Polynomial phase modelling 

A methodology for DIF estimators has been proposed based on 
polynomial phase modelling in [39] and [42]. It is assumed that the 
phase of the discrete signal may be modelled as a polynomial of order 



4>(n) = a 0 + ain + air? 4 F a p n p (54) 

The polynomial coefficients may be estimated in a number of ways. 
If the phase is unwrapped, linear least squared techniques may be ^se 



to solve for the polynomial coefficients (for this case it is assumed 
the SNR is high and the noise is additive zero-mean Gaussian). By 
generalising the result for estimating stationary tone parameters it is 
possible to derive a maximum likelihood estimator for the coefficients. 
Adaptive techniques may be incorporated into this estimator, and it 
may be implemented efficiently provided the polynomial order is small. 
Explicit details of these phase estimators as well as simulated examples 
may be found in [39] and [42]. 



6 Applications 

6.1 Automatic time- varying filtering 

The use of TFDs for non-stationary signal analysis has the advantage 
that a priori or a posteriori processing in the time-frequency plane can 
be employed to give enhanced estimates. This ability to perform time- 
frequency filtering can be applied to tasks such as signal enhancement 
and IF estimation. To perform this time-frequency filtering, a 
time- varying system transfer function, H(t, /), must be specified in 
accordance with the application. The WVD, being an entirely real 
transform, is well suited to time-varying filtering, particularly for 
monocomponent signals. The time-varying filtering operation is given 
by: 

y(t) = W- 1 [W(t,f)-H(t,f)) (55) 

where IV(£,/) is the WVD to be processed, and W~ 1 [ ] represents 
the synthesis operation. This operation determines the time domain 
analytic signal whose WVD is closest to the operand in a least mean 
squares sense [36]. If the operand is a valid WVD, the synthesis 
operator is equivalent to WVD inversion. For an in depth discussion, 
the reader is referred to chapter 17 on signal synthesis. 

For the analysis of a monocomponent signal, the time- varying filter 
would ideally capture all the signal energy and reject all the noise. In 
practice, one must settle for a tradeoff in which the significant signal 
energy is captured, while most of the noise is eliminated. A useful 
approach is to window about the IF, since it has been previously 
explained that, for monocomponent signals, the signal energy tends 
to be concentrated there. 

The various techniques proposed earlier, then, may be used for 
IF estimation as a first step in the design of the automatic filter. 
The bandwidth of the time- varying filter may be selected according to 
the application, as can the window type. For example, consider the 
problem of estimating the IF of a linear FM signal distorted by 3 dB 
white Gaussian noise. Figure 2.14 shows the IF law estimated using 
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the direct definition, while figure 2.15 shows the IF law estimated by 
using the first frequency moment of the WVD after a time-frequency 
window has been applied. In a practical situation one could start 
by obtaining a rough estimate of the IF law, then window about 
this region of time- frequency concentration, and finally obtain the 
improved IF estimate. One could also use signal synthesis to obtain 
a time domain signal from the windowed WVD, and calculate the IF 
directly from this filtered signal. 

6.2 Some specific areas of application 

In a typical seismic situation, a vibroseis signal is emitted into the 
earth. The different arrival times of reflected waveforms are used to 
reveal information as to the geological structure of the region. The 
reflected returns can also be processed to yield information as to the 
amount of dispersion present in the signal. This can be determined 
simply by comparing the IF of the received signal with that of the 
source [38] (see chapter 20 of this book on seismic application). 

In the fields of sonar and radar, the classical problem of direction- 
of- arrival estimation may be reduced to a problem of estimation of 
multiple IFs [40]. In fact, the estimation of IF is related to almost 
any problem that involves the analysis or tracking of a narrowband 
time- varying signal. This includes such natural phenomena as the 
echolocation of bats and sounds of cetacea, and more specialised 
applications such as the analysis of ECGs and helicopter blade 
rotations as presented in chapter 18 of this book by D. Forrester. 
Although the IF may not be specifically emphasised, it is very often 
the quantity of interest. 



7 Summary 

Instantaneous frequency estimation is an important problem in signal 
processing. It raises fundamental questions in signal representation, 
some of which have been addressed in this chapter. 

The original developments leading to the concept of ‘instantaneous 
frequency’ have been reviewed, and the relation of the IF to TFDs has 
been stressed. The concepts of monocomponent and multicomponent 
signals have been discussed. Techniques of IF estimation have been 
presented and some applications of IF estimation have been described. 
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